#Data:2020/09/03
#Read me:This is a sample matlab code for the Maxwell's Equations Derived
#Optimization (MEDO). 

# Contact Information:
# Email: lll_work@buaa.edu.cn

# DISCLAIMER: This sample code is provided for educational purposes only. 
# Please read the referred paper to better understand the usage and the
# inner-workings of the algorithm.
# Use at own your risk! There is no guarantee that the code is bug free.

# Please refer to the following article in your publications:
#SU Donglin, LI Lilin, YANG Shunchuan, LI Bing, CHEN Guangzhi, XU Hui, 
#A new optimization algorithm applied in electromagnetics �� Maxwell's equations 
#derived optimization (MEDO), In Journal of SCIENCE CHINA Information Sciences, 
#Volume 63, Issue 10, 2020, Pages 200301-, ISSN 1674-733X, 
#https://doi.org/10.1007/s11432-020-2927-2.
#(https://engine.scichina.com/doi/10.1007/s11432-020-2927-2)

## preset of the coefficients
#more information about the physicle meaning and the recomm range of the coefficients 
#please refer to the article above. 
def MEDO_solver():
    import numpy as np
    from utils import benchmark,derivation,solve_i
    from environment_pool import B_pool,pS_pool,g_pool,R1_pool,L2_pool,lb_pool,param_maxit_pool,popsize_pool,npar_pool,deltax_pool,ub_pool,target_couple_pool
    B = B_pool
    pS = pS_pool
    R1 = R1_pool
    g = g_pool
    L2 = L2_pool
    target_couple = target_couple_pool
    
    ## preset of the population
    popsize = popsize_pool    # population size.
    npar = npar_pool		# Dimension of the problem.
    param_maxit = param_maxit_pool		# Maximum number of iterations.
    deltax = deltax_pool #calculation accuracy of difference method

    lb = lb_pool    #Lower dimension boundary.
    ub = ub_pool  #Upper dimension boundary.

    ## the statement and initialization of the variables
    x = np.zeros((npar,popsize)) #statement of the population
    fx = np.zeros((popsize))  #statement of the objective function value
    fxmin_global = np.zeros((param_maxit))  #statement of the global optimal solution
    fxmin_f = np.zeros((param_maxit))  #statement of the optimal solution in the current iteration
    u = np.zeros((npar,popsize))  #statement of the value of the voltage sourve
    i2 = np.zeros((npar,popsize)) #statement of the current of the loop1
    i3 = np.zeros((npar,popsize)) #statement of the current of the loop2

    for i in np.arange(npar):
        x[i,:] = lb[i] + (ub[i] - lb[i]) * np.random.rand(popsize)
    #Randomize population
    x_o = 0*x #population in the last iteration, which is initialized to be 0 at the beginning.
    # print(x[:,0])
    ## calculate the objective funcion of the initialized polulation
    for i in np.arange(popsize):
        fx[i],match_request,_ = benchmark(x[:,i],npar)
    #benchmark can be replaced with your own cost function and make sure the mapping above can carry out properly.
    #print(x[:,:popsize],npar)
    xminnum = np.argmin(fx)  # finding the best solution in the initial population :
    fmin = fx[xminnum] 
    fxmin_f[0] = fmin # record the best solution in the current iteration
    xmin = x[:, xminnum] # finding the best variable in the current iteration
    glomin = fmin #the global optimal solution in the whole interations
    fxmin_global[0] = glomin #record the global optimal solution
    xminout = xmin  # record the best variable in the whole iterations


    fxlb,match_request,_ = benchmark(lb,npar) # calculate the objective function of the lower boundary.
    fmax = np.max([np.max(np.abs(fx)),np.abs(fxlb)])  # the maximum value of the objective funcion in the current iteraion
    fx=fx/fmax  # normalize the objective functions of the population.
    fxlb=fxlb/fmax # normalize the objective function of the lower boundary.
    h1 = 1 - fxlb # the variable to calculate S3
    h2 = 1 - fx[0] # the variable to calculate S3
    R2 = 0.2 # the resistant of loop1
    S3 =  0.1*np.ones((npar)) # the loop inductance of loop2 
    ## Start iterations :
    for iter in np.arange(1,param_maxit):  
        print('iter=',iter+1)
        dx = x - x_o
        x_o = x.copy()
        for i in np.arange(popsize):  
            S3 = S3 + (h1 + h2) * (x[:,0] - lb) / 2
            for j in np.arange(1,i): 
                h1 = 1 - fx[j -1]
                h2 = 1 - fx[j]
                S3 = S3 + (h1 + h2) * (x[: , j] - x[: , j-1]) / 2  
            #calculation the area of loop2: S3
            R3 = 1-fx[i] + np.sqrt((x[:,i] - lb)**2+(fx[i]-fxlb)**2) + 0.1  #calculation R3
            # the value of the voltage source is set to be the gradient of the objective function
            u[:,i] =  -derivation(x[:,i] , deltax ,npar).ravel()
            u_max=np.max(np.abs(u[:,i]))
            if(np.abs(u_max)<1e-5):
                u_max = 1e-5
            u[:,i]=u[:,i]/u_max  #normalize the voltage source
            if iter == 1: # in the first iteration, i3 is calculated in condition of static field
                Ztotal = R1 + x[:,i]-lb + R2 * R3 / (R2 + R3)
                i3[:,i] = u[:,i] / Ztotal / (1 + R3 / R2)
                i2[:,i] = R3 / R2 * i3[:,i]
                i3oo = i3[:,i] #record i3 
                i2oo = i2[:,i] #record i2
                #if(i==3):print('i3:',i3[:,3])
            else :
                # in the later iterations, i3 is calculated in condition of
                # time-varying field
                i2n,i3n = solve_i(L2,dx[:,i],1-fx[i],u[:,i],i2[:,i],i2oo,i3[:,i],i3oo,R1 + x[:,i],R2,R3,B,S3)
                i2oo = i2[:,i] #record i2
                i3oo = i3[:,i] #record i3
                i2[:,i] = i2n #update i2
                i3[:,i] = i3n #update i3
            # if(i==3):
            #     print(B,pS,g,i3[:,i],x[:,i])
            x[:,i] = x[:,i] + (i3[:,i] * B / pS - (x[:,i] - 0) * g) * (iter+1)  #update the population
        x = np.maximum(x, lb.reshape((-1,1)))
        x = np.minimum(x, ub.reshape((-1,1)))#check the population
        match_request0 = False
        for i  in np.arange(popsize):
            fx[i],match_request,_ = benchmark(x[:,i],npar)
            if(match_request):
                match_request0=True
        print('fx=',fx)
        #update the values of the objective functions
        xminnum= np.argmin(fx)  #finding the best solution in the current iteration
        fmin = fx[xminnum]
        xmin = x[:, xminnum]  #finding the best variable in the current iteration
        # if(iter==6):print(xmin)
        glomin = np.minimum(fmin, glomin) #update the global optimal solution in the whole interations
        if fmin == glomin:
            xminout = xmin # update the best variable in the whole iterations
        fxmin_global[iter] = glomin #record the global optimal solution
        fxmin_f[iter] = fmin  # record the best solution in the current iteration
        fmax=np.max([np.max(np.abs(fx)),np.abs(fxlb)])
        fx=fx/fmax  #normalization
        if(match_request0):
            break
        #print(iter+1,'x:',x[:,0])
        ## result
        # print(best_ans)
        # plot(fxmin_global) 
        # xlabel('the number of iteration')
        # ylabel('optimal solution')
        # title('convergence curve of MEDO')
        # if(iter==6):
        #     exit(-1)
    _,_,yminout = benchmark(xminout,npar)
    return xminout,fxmin_global[:iter],yminout

#C:\Users\Administrator\Desktop\Geodesic\Geodesic20.exe temp.json


#1)接收点+半径过滤网格，得到点坐标并建立目标函数
#2）接收三角形曲基函数，建立目标函数
#3）接收四边形曲基函数，建立目标函数
#4）接收NURBS曲函数，建立目标函数
#5）接收点曲线函数，建立目标函数